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Abstract. In this article we describe how two or more experimental results can be combined within the 
procedure of Feldman and Cousins, to provide combined confidence limits on the physical parameters of 
interest. We demonstrate the technique by combining the recent electron neutrino appearance results from 
T2K and MINOS. Our best fit point is sin 2 26» 13 = 0.08(0.11) and 8 = 1.1(2.0)tt; in addition we exclude 
sin 2 20i3 = at 2.7a (2.8a) for the normal (inverted) neutrino mass hierarchy. 



1. Introduction 

In order to obtain global constraints on physical parameters, it is often necessary to combine the results of 
two or more experiments with sensitivities to the same parameters. Such combinations can be performed with 
different levels of sophistication, depending on what information is available about the original measurements. 
In the absence of detailed information, relatively crude combinations, e.g. by simply summing log-likelihood 
values, can be useful, but these results will always be subject to caveats. 

In this work, we describe a method to combine results which requires somewhat detailed experimental 
information, but which focuses on obtaining correct coverage of the parameter space under study, by producing 
a combined likelihood curve from the joint data set of the experiments, and using the Fcldman-Cousins 
method [l[ to define acceptance regions for each point in the parameter space. We will describe the inputs 
that would be required from an experiment to enable its inclusion in such an analysis in ij l2.ll We also present, 
as an example of the technique, the combination of the recent electron neutrino appearance results from the 
long baseline experiments MINOS and T2K. This example is particularly pertinent, since the parameter 
013 is known to be near a physical bound; the Feldman-Cousins technique was developed to deal with such 
cases, since conventional methods for determining acceptance regions can produce incorrect coverage or null 
contours. A comparison will also be made between our method and one where fixed log- likelihood differences 
are used to produce the contours. 

1.1. The example: Constraints on #13 from MINOS and T2K 

Neutrino oscillations, a now well established phenomenon can be parametrised by the PMNS mixing 

matrix 



12l Il3j , assuming a minimal 3- flavour mixing model. All of the mixing angles in this matrix have 
now been measured to be non-zero, including the angle #13 0, Q, which governs muon to electron neutrino 
oscillations in the atmospheric (L/E) regime. The CP- violation parameter S is as yet unknown. In the 
appearance channel T2K has found an excess of electron neutrino events in a muon neutrino beam, suggesting 
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that #13 is non-zero, with a significance of 2.5a [2|. MINOS has also found an excess of events, excluding 
6*i3 = at 89% CL 0]. We will combine the results of these two experiments to get improved constraints 
on 9i3 from electron neutrino appearance measurements, as well as the CP-violation parameter S. In the 
disappearance channel, has recently been measured to be non-zero by two different reactor experiments, 
Daya Bay and RENO 

MINOS and T2K are both long baseline accelerator neutrino oscillation experiments, in which a 
predominantly beam is created by colliding protons onto a target, then magnetically focusing the resultant 
charged mesons and allowing them to decay. The neutrino beam thus produced is sampled at source by a 
near detector, and several hundred kilometers away by a far detector. Details of the individual experiments' 
setups are described elsewhere [HI, EH. Both experiments are at an atmospheric L/E, but their different 
energies and baselines (MINOS is at 735 km and has a beam peaked at 3 GeV whereas T2K is at 295 km 
with a beam peak of 0.6 GeV), along with different detector technologies and analysis techniques lead to 
different systematic errors, as well as different sensitivities to the oscillation parameters. Both experiments 
use the difference in the numbers of electron-neutrino-like events observed in their far detectors from those 
predicted to constrain #13. Both analyses are performed with the Feldman- Cousins method. 

2. Method to Combine Fits 

In order to combine the data from two or more experiments, the binned data from all experiments are 
considered together, as a single larger experiment. The Feldman- Cousins technique is then used to produce 
confidence contours based on the combined data. 

2.1. Inputs to the analysis 

In order to perform our analysis, the data from each experiment are required, along with Monte Carlo 
expectation values. Specifically, these consist of: 

(i) The number of data events in each analysis bin for each experiment. 

(ii) The expected number of events in each analysis bin, as a function of the parameters to be fitted over, 
in our example (sin 2 2#i3, 8). The number of events will vary smoothly with the oscillation parameters, 
so a set of values on a grid, which can be numerically interpolated for intermediate points, is suitable 
for this purpose. 

(iii) The correlated bin-by-bin systematic errors for each experiment, encoded in covariance matrices, as a 
function of the parameters 9. We therefore require that the calculation of a covariance matrix from 
the underlying sources of systematic uncertainty (cross-section errors, beam-line uncertainties etc.) has 
been performed by each experiment. We assume that both experimental and theoretical uncertainties 
are included in the matrix. 

Note that although we assume that the bin-to-bin correlations within each experiment will be provided, in 
general there will also be correlations between experiments, which may be important for the final result. This 
detail is discussed in § 12.61 It should be noted further that in order to be combined consistently, we require 
that the original analyses make essentially the same physics assumptions; for example the values of mixing 
parameters other than those in the fit. 

2.2. The Feldman- Cousins technique 



The Fcldman-Cousins technique is a method for generating acceptance regions within the classical framework 
for interval calculation [17[. Its distinctive feature is that when choosing which data values n to include in 



Combining Neutrino Oscillation Experiments with the Feldman- Cousins Method 



3 



the acceptance for a given 0, the candidates are ranked by the relative likelihood of observing n at 6, with 
respect to the likelihood of observing n at the best fit point for that data. The use of relative rather than 
absolute likelihoods means that data values which are unlikely for all points in the parameter space will still 
be included in the acceptance region for some 6, avoiding the problem of null contours for some data values. 
The method also automatically produces a smooth transition between one- and two-sided intervals depending 
on the observed data, with correct coverage throug 

The key procedural step in the technique is the generation of a large number of "toy" Monte Carlo 
experiments at many positions in the parameter space, in order to identify the relative likelihood limit for 
each point which will give the correct coverage. For our example, 10 4 toy experiments were generated at 
each value of considered. For each toy experiment result n, a fit is performed for #bcst using a likelihood 
function ln£. The difference, A(ln£), in the log likelihood between the best fit point ©best, and the value 
actually used to generate the toy experiment, is calculated. 

Using the A(ln£) from all toy experiments, a value A(ln£) cr i t is calculated, such that some fixed 
proportion, say 90%, of toy experiments satisfy the condition 

A(ln£) < A(ln£) crit . (1) 

The condition fl} defines an acceptance region with 90% probability for the value of 9 in question. Repeating 
the procedure for all values of 6, a surface A(ln£) cr jt(0) is calculated. 

A fit is then made to the real data using our likelihood function, obtaining a best fit point ©best, and a 
corresponding log-likelihood value (ln£)b cs t- A log-likelihood surface (ln£)(0) is also calculated on a grid in 
0. We then draw our confidence interval, using the condition ([T]), with A(ln£) = (ln£)(0) — (ln£)bost, at 
each point in 0-space to decide whether the point should be included in the contour. 



2.3. Toy experiment generation 

The toy Monte Carlo (MC) is required to give a number of events for each analysis bin, n° bs , based on 
the expected number of events n^ xp for a given value of 6, allowing for fluctuations due to systematic § 
and statistical uncertainties. The expected number of events used is the total (i.e. signal plus background) 
expectation, so that both signal and background counts fluctuate between our toy experiments. The main 
complication in the MC is to ensure that correlations between the bin values are taken into account. 

Expressing the systematic uncertainties as absolute shifts ("tweaks") xi in the expected number of events 
in each bin, we can write 

nr^n e r+ Xi . (2) 

Our method makes the common assumption that the systematic errors follow a multivariate normal 
distribution; that is, their joint pdf / sys t(£e) follows 

f syst {x) oc e-^" 1 *, (3) 

where V is the covariance matrix for the Xi. 

Since V is a symmetric, positive-definite matrix, we can use the method of Cholesky decomposition [l6| 
to find an upper-diagonal matrix L such that 

V = L T L. (4) 

| Where the number of events is very small, as with the T2K data, Fcldman-Cousins will give some over-coverage due to the 
discrete likelihood distribution generated from toy Monte Carlo events. This effect will be negligible for the combination example 
presented here. 

§ Note that this method of treating systematics is Baycsian. In a frcqucntist prescription, systematics would be included as 
extra dimensions in 6; however this is not computationally feasible when using the Fcldman-Cousins method. 
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We can then use the matrix L T as a transformation to enable us to generate the vector x from another vector 
of random variables y, such that 

x = L T y. (5) 

One can show by trivial matrix algebra that we can rewrite the factor in the exponential in (|3]) as 

x T V- 1 x = y T y, (6) 

and so we deduce that if we generate the yi independently on a normal distribution with unit standard 
deviation, then the vector x generated using ([5]) will have the desired covariance properties. From now on, 
unless otherwise stated, we assume that n^ xp has been modified for systematics using this prescription, so 
that it is now a function n^ xp (d,y) of both 8 and y. We will also use / sys t(y) to refer to the probability 
density function (pdf) of the random vector y. 

Having taken account of systematics, a number of events for each bin can be generated using Poisson 
statistics, i.e. on a pdf 

/ cxp//j \\n obs 

9pois W hs mr p (o,y)) - K jgffi ' e-^-\ (7) 

where the y{ must be generated separately for each toy experiment. 
2.4- The likelihood function In C 

When performing our fits, the value of the systematic "tweak" parameters is allowed to vary, so our likelihood 
function must include penalty terms to account for the finite uncertainty in these parameters. We define a 
likelihood function including both systematic and statistical terms, as a function of and y, using the same 
pdfs we used to generate the toy experiments: 

N 

£(n° bs ; y, 6) = [] ( ffpois « bs ; < xp (0, y))) f syst (y) 

i=0 
N 

oc II ((nT P (0,y)) n ° hS e-<^) e-^y, (8) 

i=0 

where we have dropped factors independent of (9,y). As is standard, minimisation is actually performed on 
the logarithm of the likelihood function since this is a simpler process, and is equivalent as In £ is a monotonic 
function of C. This function is given by 

N 

-21n£ = J2 (2< XP (0, V) - 2n° bs ln« x P(0, y)) + y*) . (9) 

i=0 

In the case that we do not wish to fit for systematics (e.g. for direct comparison with other analyses), we use 
the simpler function 

N 

-21n£ = 2^« xp -< bs ln< xp ), (10) 

i=0 

where the n^ xp are no longer functions of y. 
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2. 5. Performing the fit 

We perform the minimisation of (—2 ln£) using the MINUIT minimiser via the interface provided by the ROOT 
framework The derivatives c arc calculated analytically to improve performance. The fit requires 
as inputs the vector n exp , and the Cholesky decomposition L of its covariance matrix V, at every point in 9 
space. These inputs are generated on a grid in in advance and then picked out by interpolation. 

The fit is split into two nested components: a "top-level fit" over the parameters 6, and another fit 
performed over y which is called by the top-level fit to find the best (ln£) value for a given 6. sin 2 (26*13) is 
constrained to lie in the physical region [0,1]. The y are allowed to float freely, but all n° xp are constrained 
to remain above zero (actually 10 -5 ) regardless of the y. 

2.6. Correlated systematic errors 

When combining experiments with some common systematic error sources, it is necessary to calculate the 
correlation coefficients between the bins of the two experiments. In order to do this we must identify the 
common systematic uncertainties, and get from each experiment the contributions of these error sources to 
the total errors on each bin. 

Once this information is available, the cross-terms in the covariance matrix can be calculated. Assuming 
Gaussian errors, we express the number of events in each bin i as 

m = ni + cr" ncorr 2/i + CmZa ' ( n ) 

a 

where the yi and z a arc independent normally distributed random variables with (<x = 1, fi = 0), a indexes 
the correlated error sources and the a a give the contribution of the error source a to the error on bin i. 
^uncorr j g |- ne uncorrelated component of the systematic error on bin i. From (fTTj) . it is easy to show that 

V l3 = {{m -nt){ nj -nj)) = ^ CiaCja +%(C corr ) 2 - (12) 

a 

3. Our Combination Example 

In our example we combine the electron neutrino appearance results of two long baseline experiments: 
MINOS Q and T2K The inputs we take from each experiment are the expected number of electron 
neutrino events (as a function of #13, <5), in each bin, along with the covariance matrices of the systematic 
errors between the bins in each experiment. The T2K result consists of a single bin, and the MINOS result 
of fifteen bins (five of energy times three of the Library Event Matching particle identification parameter) . 

3.1. Validation 

To validate our method, we have reproduced the results of both the T2K and MINOS analyses individually. 
There are some slight differences in the ways in which the two experiments perform their analyses, which 
must be considered when validating our method, and when making a combination of their results. One 
difference is that MINOS minimise over the systematic errors in their likelihood function, whereas T2K do 
not. Since in general it may be advantageous to minimise over the systematic errors, we do minimise over 
systematics in our example fit. Another difference is that T2K show their results for sin 2 26>i3 at a fixed 
value of 2 sin 2 023 = 1, whereas MINOS show the combination 2 sin 2 6*23 sin 2 2#i3, by throwing 623, between 
its errors and then choosing #13 to keep this quantity fixed in each toy experiment. The systematic errors 
that MINOS have provided us with, however, do not include any uncertainty in the oscillation parameters 
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Figure 1. Reproduction of the T2K results: we find good agreement with the published contours 0. In 
addition, we exclude sin 2 29yj, = at 2.5o\ 

(a fixed value of 2 sin 2 023 = 1 is used). This means that our combination can be interpreted as constraining 
sin 2 20i3 at 2 sin 2 6*23 = 1. The other oscillation parameters used for calculation of the input distributions 
were: sin 2 20i 2 = 0.87, Am|i = 7.6 x 10~ 5 cV 2 and Am§ 2 = 2.3 x 10~ 3 cV 2 . For T2K, the provided inputs 
actually used Am 2 2 = 2.4 x 10 _3 cV 2 , but changing to Am§ 2 = 2.3 x 10 _3 eV 2 made negligible difference to 
the result. 

The validation results are shown in Figures Q] and [2j for T2K and MINOS respectively. In both cases, 
we find good agreement with the published contours. 

3.2. Possibly Correlated Systematic Errors 

In order to combine the errors from both experiments into a single covariance matrix, we need to consider the 
systematic errors that may be correlated between them. In the case of our example, most of the sources of 
error cited in can be assumed to be uncorrelated between experiments, due to the different experimental 
set-ups and different neutrino energies involved. However, cross section uncertainties between the single T2K 
bin (energies up to 1.25 GeV) and the MINOS low energy (1-2 GeV) bins are potentially correlated. 

To assess the necessity of evaluating this correlation, we performed a worst-case estimate of its effect, 
by assuming complete correlation between the T2K cross-section error, and the total MINOS cross-section 
plus flux error [19J. These numbers were chosen as they will overestimate the effect and are available from 
the cited publications. Running our analysis with and without these worst-case correlations included, we 
found negligible difference in the results. We therefore conclude that we can neglect correlated errors for the 
MINOS-T2K combination at this point, though this may not be true for later data sets, if cross-section errors 
assume a greater relative importance. 
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Figure 2. Reproduction of the MINOS results: we find good agreement with the published contours 

It should be noted for the future that whilst one reactor experiment could be combined with MINOS 
and T2K without needing to account for correlations (very different energies, flux, baseline etc.), if multiple 
reactor experiments were to be included some effort would be needed to calculate the correlations between 
their errors. The same would be true if another long-baseline neutrino experiment with similar peak neutrino 
energy as either MINOS or T2K were to be included. 

3.3. Results 

As previously mentioned, we take the profile likelihood, minimising over systematic errors, in our final fit. We 
also minimise over 8 when calculating the minimum value of our likelihood function, since the combination 
of the two experiments gives slight sensitivity to CP-violation. The normal and inverted neutrino mass 
hierarchies are treated separately. We find allowed regions of: 0.02(0.03) < sin 2 20 13 < 0.16(0.21) at 95% 
C.L., 0.03 (0.04) < sin 2 20 X3 < 0.15 (0.19) at 90% C.L., and 0.04 (0.05) < sin 2 20 X3 < 0.12 (0.16) at 68% C.L., 
for the normal (inverted) neutrino mass hierarchy, where we have taken the profile likelihood, minimising over 
(5, for each value of sin 2 2#i3, both for the data and during the calculation of the critical values of A(21n£) 
(<5 was thrown uniformly across [0, 27r) in the generation of the toy experiments), to give one-dimensional 
confidence intervals. Two-dimensional confidence intervals are shown in Figure [3] No values of 8 can be 
ruled out at la. The significance of the neutrino mass hierarchy preference, calculated by generating toy 
experiments about the global best fit point (which happens to be in the inverted neutrino mass hierarchy), 
and seeing what fraction of toy experiments had their global best fit point in the inverted mass hierarchy, is 
negligible. The best fit values of the oscillation parameters are sin 2 2#i3 = 0.08 (0.11) and 8 = 1.1 (2.0)7r for 
the normal (inverted) neutrino mass hierarchy. Our best fit value is compatible with the results presented 
in [ioj]. We exclude sin 2 2#i3 = at 2.7a (2.8a) for the normal (inverted) neutrino mass hierarchy. 
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For comparison, in Figure 2] we show the contours that would be obtained by selecting an allowed region 
using a fixed value of A(21n£). The regions obtained using the Feldman- Cousins approach are significantly 
narrower than those from the fixed log-likelihood contours. 

4. Conclusions 

We have demonstrated the combination of multiple experimental results with the Fcldman-Cousins method. 
Details of the inputs from experiments needed for inclusion in such fits arc outlined in § 12.11 

We would like to thank the MINOS and T2K collaborations for their help with this work. In particular 
we would like to thank: Ruth Toner and Lisa Whitehead, from MINOS, and Josh Albert, from T2K, for 
answering our many questions and helping with the validation of this work. We are also indebted to Louis 
Lyons for useful discussions. We acknowledge the support of STFC, U.K. 
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